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THE  CONDITIONAL  PROBABILITY  OF  RANDOM  HARMONIC  SETS 
Introduction 


The  approximate  probability  that  m random  variables  out  of  a set 
of  n uniformly  and  independently  distributed  on  (a,b)  are  harmonically 
related  to  within  a tolerance  of  100t7.  has  been  calculated  in  [1],  If 
accidental  harmonic  multiples  are  to  be  avoided,  the  need  for  high 
accuracy  in  frequency  measurements  is  indicated  by  these  results. 
Application  of  them  has  revealed  both  a computational  and  a conceptual 
problem  which  are  the  subjects  of  this  paper.  The  computational 
problem  is  discussed  in  connection  with  Eq.  2,  below.  A key  element 
in  its  elimination  is  a by-product  of  the  solution  of  the  conceptual 
problem:  The  derivation  in  [1]  is  valid  when  the  smallest  member  of  a 
"harmonically  related"  subset  is  a low-order  harmonic.  (It  is  compared 
there  with  the  results  of  a small-scale  simulation  for  the  cases  where 
the  smallest  member  is  the  fundamental  or  the  second  harmonic.)  If  the 
smallest  member  is  a harmonic  of  much  higher  order,  however,  asking  a 
posteriori  for  the  probabilities  of  multiples  of  the  indicated  funda- 
mental may  be  "data- snooping",  and  too  dependent  on  the  way  the  data 
fall  out.  It  is  more  appropriate  to  ask  for  the  conditional  probability 
that  all  elements  of  a subset  be  multiples  of  a fundamental,  given  that 
the  smallest  member  is  such  a multiple.  Before  addressing  this  problem, 
the  results  of  [1]  will  be  summarized  to  establish  the  notation  and  to 
disclose  the  computational  difficulty. 

Unconditional  Probability:  Summary  of  Results 

We  consider  a set  of  n random  variables  represented  by  X^,...,Xn> 
independently  and  uniformly  distributed  on  (a,b),  and  an  ordered  subset 
of  them,  Yi  £ Yo  s ...  5 Ym,  where  m=2,3,...,n.  If  the  Y^  are  all 
multiples  of  Yj/k^,  to  within  a tolerance  of  100t7.,  where  k^  is  the 
harmonic  number  of  the  smallest  member  of  the  subset,  we  say  that  they 
form  a chain  of  harmonics  of  length  m.  The  probability  of  this  is 
given  by 

Pm(ki,k2> • • « >km) 

- Pr|(Y1,...,Ym):(l-t)kiY1/k1  * Yt  * (l+t)k1Y1/k1,i-2,...,ml 


<&r(  2*> 


(i) 


where  k^  < k£  < •••  < kffl,  and  terms  of  0(t  ) have  been  dropped. 


Note:  Manuscript  submitted  May  15,  1978. 
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If  L = [bk2/a(l+t)]  is  the  largest  harmonic  permitted  by  the  range 
and  tolerance,  then  the  probability  that  a particular  subset  forms  a 
chain  of  length  m for  some  set  of  k^  is 


L-m4-3 


w 


L 


L 


k2»k2+l  k^=k2+l 


L 

L 

k =k  , +1 
m m-1 


p (k- 1 • • • »k  ) . (2) 
m l m 


The  approximate  probability  of  at  least  one  chain  of  length  m 
among  the  n random  variables  is  given  by 

^ lr 

P(n,m,k1)  = l (-l)r+1(“)  , (3) 

r-1 

where  M = (“)  . 

The  computational  problem  implicit  in  (3)  is  discussed  in  [1]  and 
is  resolved  by  a simple  recursion  formula.  A computational  problem 
not  resolved  there,  however,  is  that  of  Eq.  2.  For  a "reasonable" 
interval  (a,b),  L can  be  large,  and  the  straightforward  computation  of 
the  nested  sums  can  become  economically  prohibitive  for  m larger  than 
five  or  six.  Even  for  small  m,  repeated  calculations  can  be  costly; 
they  are  not  feasible  for  chains  of  ten  or  more. 

A partial  solution  to  this  difficulty  is  obtained  by  reformulating 
the  problem  in  terms  of  differences  between  successive  variables  rather 
than  multiples  of  a common  "fundamental".  Since  successive  differences 
can  be  the  same,  unlike  successive  multiples,  it  is  possible  to  effect 
a computational  advantage  in  an  analog  of  (2)  by  combining  terms.  This 
technique  shifts  the  computational  problem  to  longer  chains,  but  it 
does  not  solve  it.  It  is  eliminated  in  the  conditional  formulation 
which  follows. 

Conditional  Probability  of  Harmonic  Chains 

As  before,  let  Y^  s ^2  S •••  s Ym  be  an  ordered  subset  of  the  n 
random  variables  Xi,  X2,  ...  , XR  which  are  independently  and  uniformly 
distributed  on  (a,b).  The  conditional  probability  of  a chain  of 
harmonics  of  length  m,  i.e.,  the  conditional  probability  that  all 
are  multiples  of  Y./k^  to  within  a tolerance  of  100t7.,  given  that 
fQ  = Y^/k^  to  within  the  same  tolerance  is  given  by 


2 


pm(k2’“*,km  ' kl) 


Pr[|(Y1,...,Ym):(l-t)k1Y1/k1  * Y±  * (l-MOk^/k^irf m } 


0 i(Y1):(l-t)k1fQ  * Yl  * (H-t)k1fQ  } 


* ^Pr|(Y1):(l-t)k1fo  * Y1  * (l+t)^  } 


r(l+t)k1fo  m r (l+t)k.y1/k1 

-**—=  i n J 

(b-a)  (l-t)^  i=2  (l-t)k1y1/k1 

1 r>(1+t)kif0 

i dy 

(l-t)kifo 


dyi  dyl 


\0“1  m 


i»  (£E±2»  TT  k 

u v b-a.  . , i * 

i=2 


where  k^  < kj  < ***  < k^,  and  a/(l-t)  ^ k^  £ b/(l+2t) 

2 

neglecting  terms  of  0(t  ). 
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Now  let  L = [b/fQ(l+t)  ] be  the  largest  possible  harmonic.  The 
conditional  probability  that  a particular  subset  forms  a chain  of 
length  m becomes 


L-ori-2  L-m+3 


W 


/ztfor-4  v 
m*  V b-a/'  L 


k2=k1+i  k3=k2+i  v^.^1 


The  simplicity  of  (4),  as  contrasted  with  (1),  permits  an  analytic 
evaluation  of  (5),  thus  avoiding  the  computational  problem  of  (2). 
The  approximate  conditional  probability  of  at  least  one  chain  of 


length  m among  the  n random  variables  is  then  given  by  (3). 
Writing  (5)  as  a sequence  of  nested  sums,  we  have 

n ..  . , /2tfo\m_1 

W 3 {—J  s2 


where 


L-nrt-2 

V 


L-ra+3 


L-l 


*2  L 


S3*  S3  = Z 


k^k^+1 


The  innermost  sum  is 


^3*^2+^ 


/m-1 


S4’*"’Sm-l  = L 


m 


k .-k  ,+i 

m-  J.  m-2 


s = ( n k,  j [L(Lrf-l)-k  .-k  1/2 

m vi=2  1'  m-1  m-1 

Either  directly  from  (5),  or  from  (6)  and  (8),  we  have 

2tf 

W = w2  [L(L+i)-krki]  • 


(6) 


(7) 


(8) 


For  m > 2,  we  write  the  general  term  in  (7)  as 

L-i 


m-i  _ 


Sm-i+l; i=1’ ' * ’ »m“2  • 


k =*k  +1 

m-i  m-i-1 


(9) 


To  evaluate  this  sequence  of  sums  we  assume  a relationship  of  the  form 

2i+l 


m-i+1 


m-i 

(^^*)  Z Ai0(km-i)  ’ t=l,...,m-l  . 


(10) 


0=1 


For  i = m-1  this  becomes 


2m- 1 


8-1  . 


32  " Vl,0  kl  5 ”“2,3,..., 


(ID 
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from  which  (6)  can  be  evaluated  once  the  coefficients  A^^have  been 
determined.  For  i=l,  (10)  becomes 


S - ( n O [All  + A12km-1  + A13km-1 
m a=l 

Comparison  with  (8)  yields 


A^  = L(Irt-l)/2,  A^2  - A^^  = -^2» 


valid  for  all  m. 


Substituting  (10) 

in  (9)  and 

letting  k = 

k 

m-i 

we  have 

L-i 

2i+l 

,m-i-l  x 

r 

Ai»k> 

S . = 
m-i 

Ik) 

\ _ or  / 

0=1 

k=»k  . , +1 

L 

0*1 

m-i-1 

2i+l  L- 

i 

k . . 
m-i-1 

/m-i-1  \ 

v r r- 

k6  - 

V 

- 

( 11  ,k«) 

A ' 

/_  i0L  L. 

/ k0 

L- 

0=1 

0-1  k= 

■ 1 

k=l 

In  terms  of 
we  can  write 


of  the  Bernoulli  polynomials  Bn(x)  = ^ t2»  Chap*  23^ 

ite  k=o 


s 

Vi<a+1>'Vi<0) 

1 

0+1 

v 

) b 

ok  = 

k-1 

0+1 

0+1 

L. 

■y—1 

(s+l)T  , 


so  that  (13)  becomes 


, 2i+l 

i-i-1  \ V A,  q 


n 07  a#  i 


V ^ V <-  -i  ^ > 


2i+l.  0+1 


r to  i & z >r  w-z  (JiXv-i^  <i4> 


L s+l  a. 
0-1  Y-l 


We  collect  powers  of  For  6-1  we  identify  the  coefficient 

of  (k  . - 1 as  m_:L 

tn- 1 - 1 

2i+l  S+l 

A . r -i£  \ b (0Tl)[(L-i+l)Y-U;  (15) 

Ai+i,i  i~  s+l  l-  y 
0=1  Y=1 


otherwise  the  coefficient  of  (k^_^_^)  Is 


2i+l  0+1 

A . -li  (6V ^ b (S+1),  6=2,...,2i+3,  (16) 

Ai+1,5  L.  0+1  *—  ^°-l'  Y 

0-1  Y— 1 

both  for  i»l,...,m-2.  Using  these  coefficients,  (14)  becomes 


Vi  - Cv-  I WW*’1' 


which  has  the  same  form  as  (10)  with  i replaced  by  i+1.  The  values  of 
A,c  are  given  by  (12).  For  1 * i * m-2  may  be  evaluated 

recursively  using  (15)  and  (16),  with  6*1, . . . ,2m-l,  used  t 

evaluate  S2  iQ  (H)»  From  (15) 

2i+l 

Ai+l,l  * L di0Ai0  * 


where 


d.  - \ b (9+1)[(L-i-l)Y-l], 

i0  0+1  4.  Y 

Y-l 


for  i=l, . . . ,m-2  and  g=l,...,2i+l  . 


From  (16), 


Ai+1,6  " C0 6Ai0 ’ 


where 


0+1 


. Y \ (0+D 

'@6  0+1  L \b-\) y 

V=1 


(18) 


for  1=1, . . . ,m-2,  6=2... .,21+3,  e=max(l, 5-2) ,@=e, ... ,21+1,  and  C = 0 
if  0 < 6-2.  A storage  advantage  may  be  obtained  from  the  ® 

observation  that  the  may  be  stored  in  a "lower  triangular"  two- 
dimensional  array  and  the  C j may  be  stored  in  an  "upper  triangular" 
array.  Thus  they  may  share  a single  2m  by  2m- 1 array. 

The  coefficients  b^^  of  the  Bernoulli  polynomials  Bn(x)  are  given 
in  Table  23.1  of  [2]  for  0 £ n s 15.  They  may  be  obtained  in  general 
from 


.(n) 


5n_k»  k-0,...,n. 


where  Bn=Bn(0),  Bq=1,  and  [2] 

n-1 

V 


k=0 


(“)  Bk  - 0. 


\k 


The  computation  of  the  conditional  probability  of  harmonic  chains 
has  been  incorporated  into  a Fortran  IV  subroutine  called  PROHARM 
which  is  described  and  listed  in  the  Appendix.  The  required  binomial 
coefficients  are  obtained  using  subroutine  BINO,  also  listed,  as  is 
CALLPRO,  used  to  call  the  subroutines. 


APPENDIX:  Fortran  implementation 


Subroutine  PROHARM,  Figure  1,  is  a Fortran  IV  implementation  of 
the  algorithm  for  the  computation  of  the  conditional  probability  of 
harmonic  chains.  All  parameters  in  the  calling  argument  follow  the 
notation  of  this  paper.  In  principle  there  is  no  limit  to  the  length 
of  the  chain;  it  is  restricted  here  to  m ^ 20  because  of  core  limita- 
tions, with  n £ m unrestricted.  The  arrays  are  dimensioned  accordingly. 

Terms  of  order  t were  neglected  in  the  text.  Although  they  make 
no  difference  in  the  result,  it  is  not  much  more  difficult  to  calculate 
the  exact  coefficient  as 

(m-1) ! (|2_)  1 (2t )m-2 [(l+t)m  - (l-t)m] 

rather  than  the  close  approximation  ml  (2tf  //b-a))m  \ This  is  done 
on  lines  120-300  of  PROHARM.  Ll  (line  310)  is  L+l.  If  m harmonics 
will  not  "fit"  in  the  given  range,  Pr=0  is  returned  Array  AC  is  to  be 
shared  by  A^  and  Cg§,  with  the  initial  values  give,  by  (12)  stored  in 
row  22.  Throughout  the  subroutine,  the  parameter  MM  is  used  to  avoid 
repeating  a previously-completed  calculation;  only  those  coefficients 
needed  for  M > MM  are  calculated. 

The  coefficients  of  the  Bernoulli  polynomials  are  calculated 

via  lines  400-600,  and  are  stored  in  array  BB.  Array  CO  stores  the  d. 
which  are  calculated  by  lines  610-700  using  (17).  Equation  18  yields 
the  Cgj,  lines  720-840.  These  are  stored  in  the  upper  right  portion 
of  array  AC,  while  the  lower  left  of  AC  holds  the  A^  which  are 
calculated  recursively,  lines  850-950. 

Equations  11  and  6 are  implemented  in  lines  960-1010  and  1020, 
respectively,  to  obtain  Pm(ki) . Finally,  the  approximate  probability 
of  at  least  one  chain,  Eq.  3,  is  calculated  as  discussed  in  [1]  via 
lines  1030-1170.  This  calculation  usually  converges  rapidly,  and  is 
stopped  when  the  next  term  changes  the  result  by  less  than  0.0001%. 

When  the  probability  is  actually  very  close  to  unity,  however,  the 
finite  word  length  of  the  computer  may  produce  divergence  in  this 
calculation.  In  this  case,  Pr=l  is  returned.  The  validity  of  this 
procedure  has  been  established  by  observations  of  the  changes  in  the 
probability  introduced  by  small  perturbations  of  the  parameters. 

A sample  Fortran  calling  program,  CALLPRO,  is  listed  as  Fig.  2. 

This  program  retrieves  the  subroutines  PROHARM  and  BINO.  The  latter, 
used  to  calculate  the  binomial  coefficients  required  by  PROHARM,  is 
listed  as  Fig.  3.  Fig.  4 is  a run  of  CALLPRO  in  which  the  conditional 
probability  of  at  least  one  chain  of  length  m,  2 £ m £ 20,  given 
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33130 

SUBROUTINE  PPOHARM ( A,  B,  N, M,  T,  F3,  K1 , 3P.) 

33113 

DIMENSION  BBC  33, 38), AC (43, 39), C3 C IS, 37) 

33  1 23 

M 1 = M- 1 

Ij 

33  1 33 

M2*M-2 

33  1 43 

TM=  1 . 

33  1 53 

M 1 F = M 1 

33  153 

F=F1 =F3/C  3-A) 

331-73 

TMM*?D* 1 • -T 

* 

33133 

TMM=TMM*TD 

33  I 93 

TPM* TP* 1 • +T 

33233 

TOflaTPM*?*5 

33313 

IFCM.E3. 2M30T0  2 

33223 

T2=2*T 

33233 

DO  1 K=  1 , >12 

33243 

MK*M 1 - K 

33253 

TM=TM*T2 

33253 

F*F*F 1 

'?0273 

TMM»TMM*?D 

33253 

tpm=tpm*tp 

30293 

1 

Ml F*M1 F*MK 

33333 

o 

CC=M 1 F*TM*F* ( TPM-TMM ) 

33313 

LI  = INTC9/ C F3*T’5*T°)  ) + l 

33323 

IFCL1 .0E.K1 +M)G0T0  3 

33330 

PR* 3 . 

• 

33343 

RETURN 

30353 

3 

JM=2*M- 1 

33  363 

AC  (22,  1 ) = FLOAT  ( L 1 >*FL0ATCL1-1  1/2. 

33373 

IF(M . LE.MM) GOTO  6 

33353 

AC (22, 2)=AC(22, 3>  = - . 5 

33393 

IFCM.E0.2500T0  14 

30  400 

MC=2*M 1 

3341  3 

JM3=2*MM- 1 

33420 

BBC  1,  1 ) = - . 5 

33433 

13  = 3 

33  443 

IFCMM.NE.3) 10=2* (MM- 2) +3 

33453 

DO  4 1*13, MC 

03  453 

12*1-2 

33473 

Cl*l  . 

33453 

M=-  1 . 

33493 

DO  4 K* l , 1 2 

33533 

Cl=Ct*( I -K+ 1 )/K 

33513 

X=K-C1 *33( 1 , K) 

33520 

4 

BBC  1,1-1 )=v/l 

30  533 

IFCMN.EO.0)  13  = 2 

33  543 

DO  5 1=13, MC 

33553 

11=1-1 

30  563 

C 1 = 1 . 

33573 

33(1, I )= 1 . 

00583 

DO  5 K= 1,11 

33593 

C 1 =C 1 * ( I -K+ 1 1/K 

33633 

5 

BBC  I, K)  = C 1 *BB( 1, I-K) 

33613 

6 

DO  8 1=1, M2 

33623 

KM=  2* I + 1 J 

33633 

DO  8 K=  1 , KM 

33643 

lCo=K+l 

336  50 

M=3. 

30663 

FLI*FL3=FLOATCLl-I) 

33673 

DO  7 L= 1 , KP 

33683 

va’f + BB  ( KP,  L ) * ( FL I - 1 • ) 

03690 

7 

FLI=FLI*FL3 

i 
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307210 

3 

C0C I,K)=X/K? 

00  71  0 

IFCM.LE.MM)GOTO  1 1 

00720 

KM=JM-2 

00730 

MM=M 

00  740 

DO  10  J=2j JM 

00750 

J1=J-1 

00  750 

K0  = MAX3  < 1 , J-2) 

00773 

IFC J.LE. JM0)K0=JM0- 1 

00  750 

DO  10  K-K7I  j KM 

00790 

L'1=K+ 1 

00800 

X 

II 

& 

• 

00310 

DO  9 L- J 1 j LM 

00320 

CALL  3IN0CL..  JW  IG> 

00330 

9 

X=X-IG*33CLM,L) 

008  40 

1 0 

AC(J-UX)»X/LM 

00350 

1 1 

DO  13  I *2,  Ml 

00360 

KM=2*I-1 

00  3 70 

JI=XM+2 

00330 

AC  C I +2  1 * 1 ) = 0. 

00390 

DO  12  K= 1 j KM 

00900 

1 2 

ACC  1+21  ,\  )=ACC  I + 2L  1 >+CJC  I-UK)*AC(  I+23,K> 

009  10 

DO  13  J=2jJ1 

00920 

ACC  1+21 , J)=0. 

00930 

K0=MAX0C 1 , J-2) 

00940 

DO  13  K=K0 j KM 

08  9 50 

13 

AC  C 1+21 > J ) =AC  C I +21 , J ) +ACC J- 1 ,K)*ACC I +2  3/K) 

00  95  0 

1 4 

F=  1 . 

00970 

32=0  . 

00930 

FL  = FL0AT  C K 1 ) 

00990 

DO  15  J=hJM 

01000 

S2=S2+ACCM1 +21< J)*F 

01310 

1 5 

F=F*FL 

01020 

?M=CC*S2 

01030 

CALL  BIMOCMjM^  MCM ) 

01040 

KS=  1 

01050 

01  050 

IFCMCM.EQ.  1 ) RETURN 

0 1 0 70 

DO  16  K=2,MCM 

013  3 0 

K S = - 1 *KS 

01090 

S=C  CS*PM)/K)*CNCM-K+1  ) 

01130 

IFC3.LT. 1 . E-6*?R) RETURN 

01110 

IFCS.GT. 1 .E+12)G0T0  17 

01  1 20 

1 5 

?R=?R  KS*S 

01  1 30 

I F C P R • L T • 0 ) ° R=  l . 

01140 

IFCPP.GT. 1 . >PR=  1 . 

01  1 53 

RETTtrm 

01160 

1 7 

?R=  1 . 

01  1 70 

RETURN 

01  133 

END 

FIG.  1 Subroutine  PROHARM 
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30130 

PROGRAM  CALLPROC  IM°"T»  OUTPUT) 

33113 

°R I NT*  *RAMGE: 

33123 

READ#  A*  3 

301  33  1 

PRINT* * M * M * TOL*  F3*  K1 ** 

331  43 

READ*  M*  M*  T*  F3*K1 

33153 

CALL  PROHARMCA*  3*  M*  M* T*  F3*K1 * o 

33163 

ORINT  2*op 

33173 

GOTO  1 

33133  2 

FORMAT C*  5P=*( El  3 . 2/ ) 

30193 

?. ETR I EVE  C PROHARM  ) 

33233 

RETR I EVE ( B I MO  > 

33  213 

END 

FIG.  2 Program  CALLPRO 

33133 

SUBROUTINE  3IN0(H*M*K) 

33123 

L=N-M 

33  130 

MM-M 

33143 

IF(L.GT.M)MM=L 

33153 

L=N-MM 

30  163 

K=3 

33  l 70 

I F(L  • EQ  • 3 ) K=  1 

33183 

IFCL.EO. 1 )K=N 

33193 

IFCL.LE. 1 ) RETURN 

33233 

K=L 

33213 

A= 1 .*N 

33223 

D=  1 . *K 

33233 

3= A/D 

33243 

u: 

* 

CM 

II 

M 

o 

Q 

33253 

A- A- 1 . 

33263 

D=D-  1 . 

33273  1 

3=3*A/D 

33283 

X*  I NT ( 3+ . 5 ) 

33293 

return 

33333 

END 

FIG.  3 Subroutine  BINO 


’>o03°AM  CALL"’°0 


»ANjE:  A,  B ? . 35,  . 3 

N,  M,  TOL,  F3,  K1  ? 23,2.,  . 393,  .32,3 
->oa  1.33E+39 

'J,M,  TOL,  F3,  XI  ? 23,  A,  . 333,  .32,3 
otj=  9.97E-31 

M,  M,  TOL, F9,  K 1 ? 23 , 6, . 333,  . 32, g 
•'»=  g-T^E-ai 

.'J, M,  TOL,  F3,  K1  ? 23,  7,  .333,  .39,  3 
5°>  5.76E-31 

>J,**,  TOL,  F3,K1  ? 23,3,  . 333,  .32,3 
■5°*  1.33E-31 

'1,M,  TOL,  F9,  K1  ? 23,  9,  . 333,  .32,3 
1 . 97E-32 

M, M, TOL,  F9,  K1  ? 23, 1 3, .333, .32, 3 
r>D-  2.36E-33 

M,M, TOL, F3,K1  ? 23, 1 1 , . 333,  . 32, 3 
52*  1 . 57E-34 

M,  M,  TOL, F3,  K1  ? 23,  1 2, . 333, . 32, 3 
'>°=  1.35E-35 

M,M,  TOL,  F9,  KI  ? 23,  1 3, . 33 3, . 32, 3 
5®s  5.38E-37 

M, M,  TOL,  F3,K1  ? 23,  1 A,  . 333,  . 32,  A 

= 0*  1 • BSE-  3.3 

?J,M,  T0L,F3,K1  ? 23,  1 5,  .333,  .32,3 
*■5=  5.21E-13 

N, M, TOL, F3, K1  ? 23, 18, .333, .32,8 

•>9=  1.36E-11 

M,M, TOL, F3,X1  ? 23, l 7, .333, .32,3 
’»=  1.53E-13 

M, M, TOL, F3,K1  ? 23, l 3 ,. 39 3, . 32, 8 

•>»=■  2.93E-15 

N, M, TOL, F9,K1  ? 23,  19, .333,  .32,8 
5Rs  2.52E-19 

W,M, TOL, F3, K1  ? 23 , 23,  . 39 3,  . 32, 3 
00,  7.23E-2A 

FIG.  4 P(n,m,k^)  for  (a,b)  = (.05,. 8),  n=20,  t=0.003,  fo=*0.02,  k^=*8  and 

2 s m s 20. 
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f = 0.02  is  calculated.  A tolerance  of  0.003,  1^=8,  and  n=20  random 
variables  independently  and  uniformly  distributed  on  (0.05,  0.8)  are 
assumed.  With  these  parameters  and  range  of  m,  the  probability  drops 
from  unity  to  essentially  zero.  The  calculation  shown  in  Fig.  4 took 
less  than  6 CPU  seconds,  including  compilation,  on  a CDC  CYBER  174. 
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